Procesado de Señales e Imágenes Médicas

Ingeniería Biomédica

Ph.D. Pablo Eduardo Caicedo Rodríguez

2025-07-23

array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ...,
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]], shape=(200, 200), dtype=uint8)
array([[  0,   1,   2, ..., 252, 253, 255],
       [  0,   1,   2, ..., 252, 253, 255],
       [  0,   1,   2, ..., 252, 253, 255],
       ...,
       [  0,   1,   2, ..., 252, 253, 255],
       [  0,   1,   2, ..., 252, 253, 255],
       [  0,   1,   2, ..., 252, 253, 255]], shape=(200, 200), dtype=uint8)

Introduction to Wavelet Transform

Introduction

(-1.0, 6.0)
(-1.0, 5.0)

Introduction

  • It’s a mathematical tool for signal decomposition, like Fourier’s Transform.
  • Just as the Fourier transform decomposes a signal into a series of sine and cosine functions, the wavelet transform does so using a set of functions known as wavelets.
  • Wavelets are functions generated by scaling and shifting a base function known as the mother wavelet.

Introduction

Introduction

  • Morlet: Popular for time-frequency analysis in EEG and ECG.
  • Mexican Hat (Ricker): Often used in spike detection in neural signals.
  • Haar: Useful in quick decomposition of signals and feature extraction.
  • Daubechies: Frequently used in ECG signal denoising and compression.
  • Symlet: Another option for signal processing and feature extraction in EEG.
  • Coiflet: Useful for denoising and baseline correction in biomedical signals.

Introduction

  • Have a mean of zero (to capture details in the signal).
  • Be square integrable (finite energy).
  • Satisfy the admissibility condition on its Fourier transform.
  • Be oscillatory to capture frequency information.
  • (Optionally) have compact support for efficient computation and localization.

Zero Mean (Admissibility Condition)

The function must have an average value of zero. Mathematically, this is expressed as:

\[\int_{-\infty}^{\infty} \psi(t) \, dt = 0\]

This condition ensures that the wavelet can detect changes or “details” in the signal rather than its average or constant components.

Square Integrability

The function \(\psi(t)\) must be square integrable, meaning it has finite energy:

\[\int_{-\infty}^{\infty} |\psi(t)|^2 \, dt < \infty\]

This requirement ensures that the wavelet’s energy is finite, making it possible to localize the function in both time and frequency domains. Functions that satisfy this belong to the \(L^2(\rm I\!R)\) space, which is the space of all functions with finite energy.

Admissibility Constant

The wavelet’s Fourier transform, \(\hat{\psi}(\omega)\), should satisfy the admissibility condition:

\[C_\psi = \int_{-\infty}^{\infty} \frac{|\hat{\psi}(\omega)|^2}{|\omega|} \, d\omega < \infty\]

where \(\hat{\psi}(\omega)\) is the Fourier transform of \(\psi(t)\), and \(\omega\) represents angular frequency. This condition implies that \(\hat{\psi}(\omega)\) must approach zero as \(\omega \rightarrow 0\) meaning the wavelet has no component at zero frequency (or DC component). This condition is crucial for ensuring that the wavelet transform is invertible.

Oscillatory Nature

A mother wavelet should generally be oscillatory or “wavelike” (hence the term “wavelet”). This oscillatory behavior allows the wavelet to capture variations in the signal. For example, wavelets like the Morlet wavelet resemble decaying sinusoids. This oscillatory nature helps the wavelet capture both high-frequency and low-frequency components effectively.

Compact Support

Although not strictly necessary, compact support is often a desirable property. Compact support means that the function is non-zero only over a finite interval, making it well-localized in time. This allows for efficient computation and good localization in the time domain. For example, the Haar wavelet has compact support, while others, like the Morlet wavelet, do not have strict compact support but still decay rapidly.

A wavelet transformation

(0.0, 1.0)

A wavelet transformation

Mathematical Expressions

The continuous wavelet transform of a signal \(f(t)\) is defined as:

\[W_f(a, b) = \int_{-\infty}^{\infty} f(t) \, \psi^*\left(\frac{t - b}{a}\right) \, dt\]

where:

  • \(f(t)\) is the input signal,
  • \(\psi\) is the mother wavelet,
  • \(a\) is the scale parameter (controls the width of the wavelet),
  • \(b\) is the translation parameter (controls the position of the wavelet),
  • \(\psi^*\) denotes the complex conjugate of the mother wavelet. ## Mathematical Expressions

The continuous wavelet transform of a signal \(f(t)\) is defined as:

\[W_f(a, b) = \int_{-\infty}^{\infty} f(t) \, \psi^*\left(\frac{t - b}{a}\right) \, dt\]

where:

  • \(f(t)\) is the input signal,
  • \(\psi\) is the mother wavelet,
  • \(a\) is the scale parameter (controls the width of the wavelet),
  • \(b\) is the translation parameter (controls the position of the wavelet),
  • \(\psi^*\) denotes the complex conjugate of the mother wavelet.

Mathematical Expressions

The inverse continuous wavelet transform is given by:

\[f(t) = \frac{1}{C_{\psi}} \int_{0}^{\infty} \int_{-\infty}^{\infty} W_f(a, b) \, \psi\left(\frac{t - b}{a}\right) \frac{db \, da}{a^2}\]

where:

where \(C_{\psi}\) is a normalization constant, defined as:

\[C_{\psi} = \int_{0}^{\infty} \frac{|\hat{\psi}(\omega)|^2}{\omega} \, d\omega\]

and \(\hat{\psi}(\omega)\) is the Fourier transform of the wavelet \(\psi(t)\).

Mathematical Expressions

The discrete wavelet transform decomposes the signal at discrete levels of scale. For a signal \(x[n]\), the wavelet decomposition is defined as:

\[c_{j, k} = \sum_{n} x[n] \, \psi_{j, k}(n)\]

where:

  • \(\psi_{j, k}(n)= \frac{1}{\sqrt{2}}\psi\left(\frac{n-2^{j}k}{2^{j}}\right)\) represents the scaled and translated versions of the mother wavelet \(\psi\)
  • \(j\) is the scale index, and \(k\) is the translation index.

The decomposition typically consists of approximation and detail coefficients at each scale:

Approximation coefficients \(a_j\): capture the low-frequency components. Detail coefficients \(d_j\) capture the high-frequency components.

Mathematical Expressions

The inverse discrete wavelet transform reconstructs the original signal from its approximation and detail coefficients: \[x[n] = \sum_{j} \sum_{k} c_{j, k} \, \psi_{j, k}(n)\]

This reconstruction process involves upsampling and filtering of the approximation and detail coefficients at each scale.

Using CWT

  • Purpose: The CWT is used when you need a highly detailed, continuous analysis of a signal over all possible scales and positions.
  • Output: CWT gives you a “heatmap” of wavelet coefficients, showing which frequencies (or scales) are present in the signal at each point in time. This allows for a continuous representation.
  • Applications: CWT is useful for analyzing signals where you want to see the evolution of frequencies over time, such as:
    • Detecting subtle changes in frequencies over time (like brainwave analysis in EEG).
    • Signals with non-repeating, transient features (like spikes in biomedical signals, e.g., ECG).
  • Trade-Off: CWT is more computationally expensive because it analyzes all scales and translations continuously. It gives lots of data but is slower and requires more memory.

Use CWT when:

  • You need a detailed, continuous representation.
  • You want to detect subtle or fast-changing features across time.
  • You’re okay with higher computational costs to get very fine-grained analysis.

Using DWT

  • Purpose: The DWT is used when you want a compact, efficient representation of a signal, usually for compression or feature extraction. It analyzes only specific scales (powers of two), not continuously.
  • Output: DWT gives you a set of coefficients at each level (or scale), capturing information at that specific scale. It’s efficient and uses fewer data points.
  • Applications: DWT is ideal when you want to reduce the size of data or focus on a smaller set of frequencies, such as:
    • Image and audio compression (like JPEG 2000 or MP3 formats).
    • Feature extraction for machine learning (e.g., identifying specific patterns).
    • De-noising signals by discarding certain scales that contain noise.
  • Trade-Off: DWT is computationally cheaper but less detailed than CWT. It doesn’t give a continuous heatmap but rather a discrete set of scales.

Use DWT when:

  • You need a compact and efficient representation.
  • You’re focused on data compression, de-noising, or feature extraction.
  • You want faster computations with less data storage requirements.

Multiresolution Analysis (MRA) of ECG

Introduction: Analyzing Complex Signals

  • Signals like ECG are often non-stationary.
  • Their frequency content and characteristics change over time.
  • Traditional methods (e.g., Fourier Transform) analyze the signal globally, losing temporal information.
  • Need a method to analyze signals at different time scales and frequency bands simultaneously.

What is Multiresolution Analysis (MRA)?

  • A mathematical framework to decompose a signal into components at different levels of resolution or detail.
  • Think of it like looking at a picture:
    • Coarse view (low resolution): See the overall structure.
    • Fine view (high resolution): See small details.
  • MRA decomposes a signal into a series of approximations and details, capturing information across various scales.
  • Often implemented using Wavelet Transforms (specifically the Discrete Wavelet Transform - DWT).

MRA: The Nested Subspace Concept

  • MRA is built upon a sequence of nested subspaces \(V_j\) in \(L^2(\mathbb{R})\).
  • These subspaces represent approximations of the signal at different resolutions.
  • Key Properties:
    • Nesting: \(\dots \subset V_2 \subset V_1 \subset V_0 \subset V_{-1} \subset V_{-2} \subset \dots\)
    • Density: The union of \(V_j\) is dense in \(L^2(\mathbb{R})\).
    • Separation: The intersection of \(V_j\) is \(\{0\}\).
    • Scaling: \(f(t) \in V_j \iff f(2t) \in V_{j-1}\)

The Scaling Function (\(\phi\))

  • Each subspace \(V_j\) is generated by scaled and translated versions of a single function called the scaling function, \(\phi(t)\).
  • \(\{\phi(t-k) : k \in \mathbb{Z}\}\) forms an orthonormal basis for \(V_0\).
  • \(\{\sqrt{2^j}\phi(2^j t-k) : k \in \mathbb{Z}\}\) forms an orthonormal basis for \(V_j\).
  • The scaling function captures the “smooth” or low-frequency components – the approximation of the signal at a given scale.

The Wavelet Function (\(\psi\))

  • The difference in information between two successive approximation spaces (\(V_j\) and \(V_{j-1}\)) is captured by the detail spaces, \(W_j\).
  • \(V_{j-1} = V_j \oplus W_j\), where \(W_j\) is the orthogonal complement of \(V_j\) in \(V_{j-1}\).
  • There exists a function \(\psi(t) \in W_0\), the wavelet function (or mother wavelet).
  • \(\{\psi(t-k) : k \in \mathbb{Z}\}\) forms an orthonormal basis for \(W_0\).
  • \(\{\sqrt{2^j}\psi(2^j t-k) : k \in \mathbb{Z}\}\) forms an orthonormal basis for \(W_j\).
  • The wavelet function captures the “details” or high-frequency components.

Linking \(\phi\) and \(\psi\): Two-Scale Equations

  • The scaling and wavelet functions are related through coefficients \(h_k\) and \(g_k\).
  • \(\phi(t) = \sqrt{2} \sum_k h_k \phi(2t - k)\)
  • \(\psi(t) = \sqrt{2} \sum_k g_k \phi(2t - k)\)
  • These equations show how functions at one scale (left side) are constructed from functions at a finer scale (right side).
  • The coefficients \(h_k\) and \(g_k\) are critical – they define the specific wavelet and are the impulse responses of filters used in discrete implementations.

From Math to Practice: Filters

  • The sequences \(h_k\) and \(g_k\) correspond to digital filters:
    • \(H = \{h_k\}\) is a low-pass filter.
    • \(G = \{g_k\}\) is a high-pass filter.
  • These filters are designed such that \(G\) is derived from \(H\) (for orthonormal wavelets, \(g_k = (-1)^k h_{N-1-k}\)).
  • Applying the low-pass filter and downsampling corresponds to projecting the signal onto \(V_j\).
  • Applying the high-pass filter and downsampling corresponds to projecting the signal onto \(W_j\).

Discrete MRA: Mallat’s Algorithm

  • The discrete implementation of MRA is efficiently computed using Mallat’s Algorithm (or the Fast Wavelet Transform - FWT).
  • It’s a pyramidal algorithm based on filter banks and downsampling/upsampling.
  • Decomposes the discrete signal into approximation (\(A\)) and detail (\(D\)) coefficients at successive levels.

Mallat’s Algorithm: Decomposition

Starting with signal \(x[n]\) (considered \(A_0\)):

  1. Convolve \(A_j\) with low-pass filter \(H\) and high-pass filter \(G\).
  2. Downsample outputs by 2 (\(\downarrow 2\)).
    • Output of \(H \to \downarrow 2\) gives approximation coefficients \(A_{j+1}\).
    • Output of \(G \to \downarrow 2\) gives detail coefficients \(D_{j+1}\).
  3. Repeat the process on the approximation coefficients \(A_{j+1}\) for the next level.

\[ A_{j+1} = (A_j * H) \downarrow 2 \] \[ D_{j+1} = (A_j * G) \downarrow 2 \]

Mallat’s Algorithm: Reconstruction

Reconstructing signal \(A_j\) from \(A_{j+1}\) and \(D_{j+1}\):

  1. Upsample \(A_{j+1}\) and \(D_{j+1}\) by 2 (\(\uparrow 2\), insert zeros).
  2. Convolve upsampled coefficients with reconstruction filters \(H'\) and \(G'\).
  3. Add the results to get \(A_j\).
  4. Repeat until the original signal \(A_0\) is reconstructed.

\[ A_j = (A_{j+1} \uparrow 2 * H') + (D_{j+1} \uparrow 2 * G') \] (For orthonormal wavelets, \(H'\) and \(G'\) are time-reversed \(H\) and \(G\))

MRA Architecture

Taken from Stackexchange

Why MRA for ECG?

  • ECG signals are non-stationary and have features at different scales:
    • QRS complex: Sharp, high frequency, short duration.
    • P and T waves: Broader, lower frequency, longer duration.
    • Baseline wander: Very low frequency.
    • Muscle noise: High frequency.
  • MRA naturally separates these components into different frequency bands (detail levels), making analysis easier.

Applications of MRA in ECG: Denoising

  • Noise (powerline, muscle, baseline wander) occupies different frequency bands.
  • Decompose noisy ECG using MRA.
  • Noise components are isolated in specific detail coefficients:
    • High-frequency noise in fine detail levels.
    • Baseline wander in coarse approximation/detail levels.
  • Apply thresholding to the noise-dominated coefficients (e.g., set small values to zero).
  • Reconstruct the signal from the modified coefficients to get a denoised ECG.

Applications of MRA in ECG: Feature Extraction

  • Different ECG waves are highlighted in different MRA levels:
    • QRS complex: Often prominent in mid-to-high frequency detail coefficients.
    • P and T waves: Appear in lower frequency detail coefficients or approximation coefficients.
  • Analyze coefficients at relevant scales:
    • Detect QRS peaks by finding local maxima in specific detail levels.
    • Delineate P and T waves by analyzing the shape and zero-crossings of coefficients at coarser scales.
  • Extract clinically relevant features (intervals, amplitudes) from the detected waves.

Other Applications in ECG

  • Compression: MRA provides a sparse representation; many coefficients are small and can be discarded or quantized, reducing data size.
  • Abnormality Detection and Classification: Features derived from the distribution or patterns of coefficients across different MRA levels can be used as input for machine learning algorithms to classify arrhythmias or other heart conditions.

Practical Implementation Steps

  1. Acquisition & Preprocessing: Get digitized ECG data. Basic filtering might precede MRA if necessary.
  2. Choose Wavelet: Select a wavelet family (e.g., Daubechies ‘dbN’, Symlets ‘symN’). Consider similarity of wavelet shape to ECG features (e.g., QRS) and vanishing moments.
  3. Determine Decomposition Levels: Choose the number of levels (\(J\)) based on sampling frequency (\(F_s\)) and the frequency bands of interest. Detail level \(j\) covers frequencies approx. \([F_s/2^{j+1}, F_s/2^j]\).
  4. Perform Decomposition: Apply Mallat’s algorithm (DWT) to get approximation and detail coefficients.
  5. Analyze/Process Coefficients: Apply denoising (thresholding) or feature extraction techniques to the relevant coefficients.
  6. Reconstruct (Optional): Reconstruct the signal if a filtered or modified ECG is needed.

Implementation: Wavelet Choice & Levels

  • Wavelet Family: Daubechies (db) and Symlets (sym) are common for ECG. ‘db4’ is often cited for its similarity to the QRS complex.
  • Number of Levels (\(J\)):
    • Higher levels give coarser approximations and lower frequency details.
    • Choose \(J\) so that key ECG features fall into specific, separable detail levels.
    • E.g., for a 360 Hz ECG, QRS (10-50 Hz) might be in D3-D5, P/T (1-10 Hz) in D5-D7, baseline wander (<1 Hz) in A7 or D8+.

Conclusion

  • Multiresolution Analysis is a powerful technique for analyzing non-stationary signals like ECG.
  • It decomposes the signal into different frequency bands/scales using scaling and wavelet functions.
  • Implemented efficiently via Mallat’s algorithm (DWT).
  • Enables effective ECG denoising by isolating noise at specific levels.
  • Facilitates robust feature extraction (P, QRS, T) by highlighting them in different detail coefficients.
  • A fundamental tool in modern automated ECG analysis systems.

ECG Analysis

# Elegir wavelet y nivel de descomposición
wavelet = "db4"
max_level = 3

# Descomposición multiresolución
coeffs = pywt.wavedec(ecg001, wavelet, level=max_level)
# coeffs[0]: coeficientes de aproximación al nivel 3 (cA3)
# coeffs[1:], coeficientes de detalle en niveles 3, 2, 1 (cD3, cD2, cD1)

# Reconstruir aproximación al nivel máximo
arr_approx = [coeffs[0]] + [np.zeros_like(c) for c in coeffs[1:]]
approx = pywt.waverec(arr_approx, wavelet)

ECG Analysis


# Reconstruir detalles por cada nivel
details = []
for i in range(1, max_level + 1):
    arr_detail = [np.zeros_like(c) for c in coeffs]
    arr_detail[i] = coeffs[i]
    detail = pywt.waverec(arr_detail, wavelet)
    details.append(detail)

ECG Analysis

ECG Analysis

ECG Analysis